Brownian motion with dry friction: Fokker-Planck approach 
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We solve a Langevin equation, first studied by de Gennes, in which there is a solid-solid or dry 
friction force acting on a Brownian particle in addition to the viscous friction usually considered in 
the study of Brownian motion. We obtain both the time-dependent propagator of this equation and 
the velocity correlation function by solving the associated time-dependent Fokker-Planck equation. 
Exact results are found for the case where only dry friction acts on the particle. For the case where 
both dry and viscous friction forces are present, series representations of the propagator and corre- 
lation function are obtained in terms of parabolic cylinder functions. Similar series representations 
are also obtained for the case where an external constant force is added to the Langevin equation. 

PACS numbers: 05.40.-a, 02.50.-r, 05.70.-a 

I. INTRODUCTION 

In one of his last papers, P.-G. de Gennes [1 proposed to study Brownian motion under the influence of a solid- 
solid or dry friction force 1 in addition to the viscous friction force commonly studied since the works of Einstein and 
Langevin [3J. The main property of dry friction, as is known from common experience, is that a certain threshold 
force Ap > has to be applied to a solid object resting on a solid surface in order to move the object [2]. If the 
applied force F is smaller than Ap, then the object stays at rest (stick state), whereas if F > Ap, then the object 
starts moving (slip state). 

The simplest way to model this behaviour mathematically is to consider a velocity-dependent force having the form 
—Apa(v), where <r(v) stands for the sign of v with the convention tr(0) = 0. By incorporating this force into Newton's 
equation, 

mv = —Apa(v) + F, (1) 

it is indeed easy to see that v = is the unique attracting state when \F\ < Ap, and that v = ceases to be attracting 
when \F\ > Ap. An object of mass m initially at rest, v(0) = 0, will therefore remain at rest unless an external force 
F greater in magnitude than the threshold or contact force Ap- is applied. 

This simple description of dry friction is commonly used in engineering and physics to study macroscopic systems 
in which solid-solid friction plays a significant role @H7]. It is used, for example, to study stick-slip motion [8HTT]. 
vibrated granular media |12l 113] , as well as the ratchet motion of solid objects moving on a surface which is vibrated 
horizontally by a time-dependent force [14] (see also [E]). In this context, de Gennes PQ proposed to study the case 
where the surface vibrates randomly by considering F in ([!]) to be a Gaussian white noise. Under this assumption, he 
derived an approximate expression for the velocity correlation function (v(t)v(0)) which generalises the well-known, 
exponentially-decaying velocity correlation function of the linear Langevin equation with viscous friction [3j. He also 
obtained the stationary velocity distribution, previously reported in (TSJ US] . 

In this paper, we extend de Gennes's work by calculating the time-dependent propagator p(v, t\vo, 0) of the Langevin 
equation 

mi = -av - Ap<j(v) +F + £(*), (2) 

which includes a viscous-like friction force with drag a, a dry friction force with threshold force Ap, a constant 
external force F, and a Gaussian white noise £(t) characterised by 

«(*))= o, mm) = ™ 2 rs(t), (3) 

with m 2 T the noise power. In [El [18], the propagator of ^ was obtained using path integral techniques in the small- 
noise limit, r — > 0. Here we obtain the propagator without any approximation by directly solving the time-dependent 
Fokker-Planck equation associated with ([2]), which is explicitly given by 

dp{v,t\v o ,0) d . r d 2 p(v,t\v ,0) 

Wt = dv^ v + (u) ~ °M«.*l«o,o) + 2 > ( 4 ) 



Dry friction is also called Coulomb friction [2]. 
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where 7 = a/m, A = Ap/m, and a — F/m denote the force parameters rescaled by the mass m. To approach the 
solution of the Fokker-Planck equation, we first consider in section [n] the case where only dry friction is present, i.e., 
7 = F = 0, in order to establish the notation and to restate the results of de Gennes pQ. The effects of the viscous 
force is then studied in section |III| In section |IV| we finally consider the case where all forces are present in order to 
study the stick-slip transition described above, but now in the presence of noise. We shall see there that, although 
there is strictly speaking no sticking in the presence of Gaussian white noise pQ, the Langevin dynamics with dry 
friction does display two phases that resemble the stick and slip states of the noiseless dynamics. 

Apart from extending de Gennes's results to include viscous damping and external forcing, the model studied here 
is interesting in part because it provides an example of piecewise-smooth dynamical system perturbed by noise [H] 
which can be solved exactly by analytical means. In the case of pure dry friction, we are indeed able to obtain closed- 
form expressions for the time-dependent propagator and the velocity correlation function, while for the dry-viscous 
and full-force cases, these functions are written in terms of spectral sums involving the so-called parabolic cylinder 
function. These exact results are valid for any noise power, and so also extend those obtained in [T71 118j . 

In a sense, it should not come as a surprise that the Fokker-Planck equation Q can be solve exactly, since it 
is one-dimensional and piecewise linear. However, the solution is non-trivial and provides some physical insights 
into stick-slip transitions perturbed by noise. In particular, we are able to identify from it a number of time scales, 
expressed in terms of the noise power and force coefficients, which could be valuable inputs for experiments on driven 
systems with dry friction and driven droplet dynamics, such as those performed recently by Chaudhury et al. |20H23) . 
Our results on the spectrum of the Fokker-Planck operator of Q also enable us to state conditions under which the 
viscous force can be neglected with respect to the dry friction, a regime referred to by de Gennes pQ as the partly stuck 
regime, or, conversely, under which dry friction can be neglected with respect to viscous friction (viscous regime). 
These conditions might be useful for determining whether dry friction or friction-like forces, such as liquid hysteresis 
[20] [21], should be considered when modeling noise-driven systems. 



II. DRY FRICTION ONLY 



We recall in this section the solution of the Fokker-Planck equation Q in the simple case where there is no viscous 
damping and no external force, i.e., 7 = 0, a = 0, A>0 and T > 0, for the purpose of introducing the general 
notation used throughout the paper and recalling the results of de Gennes PQ. 

The solution proceeds in the standard way (see, e.g., [24] and [25] for a rigorous but accessible account). Using the 
appropriate nondimensional variables 

2A 2A 2 

x=yrv, T =—t, (5) 

we first recast the Fokker-Planck equation in the form 

dp(x,r\x',0) d<S>'{x)p(x,r\x',0) d 2 p(x, t\x', 0) 



Or dx dx 2 



(6) 



where x' stands for the initial condition, and $(x) — \x\ is the nondimensional potential associated with the dry 
friction force. The normalised stationary solution of |6]) is given in terms of this potential by 

e -*( a ) 

P*{x) = — — , (7) 

where Z* is the normalisation factor. 

To find the time-dependent state p(x, t\x', 0), we next solve the eigenvalue equation 

- Au A (x) = (&{x)u A (x))' + ul(x), (8) 

for A > 0. The eigenfunctions v A of the adjoint problem are related to the eigenfunctions u A via the stationary 
distribution 

u A (x) = v A (x)e-*M. (9) 
The two sets of eigenfunctions form a bi-orthogonal system, so that the propagator can be formally written as 

p(x, t\x', 0) = J2 e- AT u A (x)v A (x')/Z A (10) 
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FIG. 1: (Color online) Propagator p(x, t\x' , 0) of the Fokker-Planck equation for the pure dry friction case with initial condition 
x' — 2. The dashed line shows the stationary solution. 



with normalisation 



v A (x) u a (x) dx. 



(11) 



Since we are concerned with a piecewise smooth potential whose derivative, representing the force, shows a finite 
jump at the origin, we have to solve Q on the positive and negative domains separately, with standard, decaying 
boundary conditions at infinity, and then match at the origin the two solutions obtained by imposing continuity of 
the eigenfunction and continuity of the probability current 



This yields the conditions 



j(x) = <S>'(x)u A {x)+u A (x). 



ua(O-) - u A (0+) 

$'(o-)u A (o-) + <(o-) = $'(o+K(o+) + *4(o+). 



(12) 



(13) 



Additionally, we impose the initial condition p(x, 0\x' , 0) = 6(x — x'). 

The complete analytical solution of ^ with the above boundary conditions can be easily written down; see, e.g., 
[26) for a classical account. To keep the presentation self-contained, we summarise in Appendix [A] the major steps 
leading to the final solution given by 



p(x,r\x',0) 



-t/4 



2y/TTT 



\x\~\x'\)/2-(x-x') 2 /(iT) e 1 



+ 



1 + erf 



(M + W\) 

2^ 



(14) 



where erf(ar) = 2 L exp(— z 2 )dz/y/ir denotes the probability integral. 

Figure [I] shows the behaviour of this solution as r increases. At a qualitative level, the behaviour of p(x, t\x', 0) 
shows two stages: a short time relaxation of the initial distribution towards the origin with the build-up of a corner 
or cusp at the origin, followed, on a longer time scale, by a build-up of the exponential tails of the stationary state 
p*{x). The first stage can be understood superficially by noting that the deterministic (zero- noise) dynamics of the 
model admits the solution v(t) — vq — At with v$ > 0, and thus shows a finite-time relaxation towards the fixed 
point v = after a time t c — vq/A. In the nondimensional variables of ([5]), this time-scale reads r c = x' and is the 
time-scale within which the propagator develops its unimodal, exponential-type shape (see figure [lj. 

The effect of the noise on the finite-time decay of v(t) can be seen by computing the expectation of x(r), i.e., the 
first moment of ( 14 ) : 



(x(r)) = a(x') 



_ erfc(-(|^|-r)/(2V?)) + , + r) exp(| a ;1)erfc((|^|+r)/(2^)) 



(15) 



Here erfc(x) = 1 — erf (a;) denotes the complementary error function. If one considers cases with large initial conditions 
(\x'\ ^> 1), i.e., cases where the noise can be considered as a small perturbation, then the initial stage r < \x'\ is 
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dominated by the first contribution to ( |T5| , and the mean value follows the deterministic path x' — ra{x'), apart from 
exponentially small corrections. Beyond this time scale, i.e., r > |sc'|, both contributions are of the same size, and the 
mean value follows an exponential relaxation towards its stationary value X* = 0, dominated by diffusive behaviour. 
To complement this analysis, we can compute the correlation function, defined by 

(x(r)x(O)) = / dx dx'xx'p(x,T\x',0)p*(x'). (16) 



Because of the spectral gap observed in the spectrum of the Fokker-Planck operator (see Appendix [A}, we know that 
(x(t)x(O)) must decay exponentially for large times; however, because of the finite-time decay of the dry friction 
dynamics, we expect to see corrections to this exponential decay for short times. This can be verified explicitly by 
computing the two integrals appearing in the expression of (x(t)x(O)) and by noting that the symmetric part of the 



propagator (14) does not contribute in these integrals. The result is 



(x(t)x(O)) = —= / dx I dyxye- ( - x+v ^ 2 e-^ +v ^^ AT hmh(xy/(2T)) 



o 



e~ r/4 I / V7TT 



(, { (^p eT/4erfc (\^/ 2 ) - i) P + 6r2 - 12r + 24 ) + 2t2 + 24 | 

where, in the last step, the standard asymptotic property of the probability integral was used. As usual, we can define 
a correlation time r corr from this result by putting the dominant exponential term in the form e~ T / Tcorr to obtain 

Tcorr = 4 Or 

- 2m2r 

^corr — T~2 ( "J 

F 



in the original variables given in (15]). This correlation time and the last line of (17) recover the results of de Gennes 
PQ up to some constants. 2 

III. DRY AND VISCOUS FRICTION 

For the case with viscous friction, 7 > 0, but without external force, a = 0, the appropriate nondimensional units 
are given by 

27 x 1/2 



x = [ jr ) v, t = jt. (19) 



In terms of these variables, the Fokker-Planck equation Q takes the form ^ with the potential 



2 



1/2 



$(x) = (|x|+<5) 2 /2, S = A [^) ■ ( 2 °) 

The parameter S measures the strength of the dry friction relative to the viscous damping and is therefore positive. 
However, a negative 5 is possible if we have in mind to study Kramer's classical problem of a Brownian particle 
diffusing in a bistable or double-well potential. We will mention this case later in the paper. 



Since the potential (20) is piecewise parabolic, the eigenvalue problem of the Fokker-Planck operator can still be 
solved analytically. However, to the best of our knowledge, a closed analytical expression of the propagator is not 
known. Various details can be found in the literature, in particular, in the context of the exit-time problem for the 
Ornstein-Uhlenbeck process [27] and the constrained harmonic quantum oscillator 28, 29 . As a coherent presentation 
cannot easily be found, we provide some details in Appendix |B| 



2 There is a typo in equation (37) of pQ: tA should be t/\. There also seems to be a missing factor 2 in some of de Gennes's results. 
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FIG. 2: (Color online) (a) Parabolic cylinder function Da{S) as a function of the index A for 5 = (purple) and 5 = 0.5 
(blue). The zeros yield the eigenvalues a£°' which, for 5 = 0, are the odd positive integers, (b) Parabolic cylinder function for 
5 = 5 normalised by the asymptotic growth Dasymp = 2 A/2 P((1 + A)/2) for large A; see [30]. The resulting expression remains 
bounded. 
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FIG. 3: (Color online) Eigenvalues of the dry friction and viscous friction problem as a function of 8; see (21 1 and (22 1. The 
red line corresponds to Ao = 0, while the solid and dashed lines correspond, respectively, to the odd and even eigenvalues for 
n = 1, 2, 3, 4, and 5 (see also [2S] H2 [31]). 



Because of the quadratic shape of the potential (20), we expect the spectrum to be pure point. The eigenfunctions 

associated with the odd 



are either even or odd because of the symmetry of the potential. The eigenvalues A 
eigenfunctions are determined for n > 1 by the characteristic equation 



(o) 



0. 



(21) 



where D\(z) denotes the parabolic cylinder function [30] (see figure |2h. 4 Apart from the stationary solution with 
eigenvalue Aq = 0, the other eigenvalues associated with the even eigenfunctions are simply given by the relation 



A( e ) = A (o) + 1. 



(22) 



The case without dry friction, i.e., 5 — 0, results in a quadratic potential associated with the Ornstein-Uhlenbeck 
process. Even and odd eigenvalues are then given by positive even and odd integers, as the characteristic equation 
(21| reduces to an equation involving Hermite polynomials. For negative values of 5, the potential becomes bistable 
and the lowest non-trivial eigenvalue approaches zero, reflecting the behaviour of Kramer's escape rate |24j . For large 



s Elliot's Theorem |25| does not apply here, so a rigorous proof of the pure point property cannot easily be given. 

4 As expected, all the eigenvalues are positive, since the parabolic cylinder function is strictly positive for negative values of the index A. 
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FIG. 4: (Color online) Eigenvalues for the dry an d viscous friction case, (a) Illustration of the asymptotic relation (lines) 
satisfied by A^ ' (symbols) for large index n; see {gg} , (b) Ai o) (see also fi gure|2J) as a function of 8 together with the asymptotic 
expressions (23 1 (gray dashed line) and (24 1 (black dashed line). A crossover between the small 8 (large n) and large 8 (small 
n) asymptotics is clearly visible, (c) Log-log plot of the difference between the eigenvalue a£°' and the asymptotic expression 
(23 1, referred to as Ajf ymp in the plot, as a function of n. The slope of the observed straig ht l ines is approximately —1/2. (d) 
Log-log plot of the distance between the eigenvalue An°' and the asymptotic expression of (241, referred to as A^ symp (i5) in the 
plot, as a function of 8. The slope of the observed straight lines is approximately —2. 



positive values of 5, the spectrum develops a gap between the two lowest eigenvalues, Aq and , and the rest of the 
spectrum. These properties of the characteristic equation and of the spectrum are illustrated in figures [2] and [3] 

Given the closed analytic expression of the characteristic equation, we can use standard asymptotic expansions 
for the parabolic cylinder function to obtain estimates for the eigenvalues |30) . In this respect, the expansion of the 
oscillating part of the parabolic cylinder function, e.g., in terms of Airy functions, is quite useful, as it enables us to 
obtain expressions for the eigenvalues for large index n and for large values of the parameter 8. For the upper part of 
the spectrum, we explicitly obtain 



28 



A^ ~ (2n - 1) + — V2n- 1 + 2 



n > 1. 



(23) 



This asymptotic relation is illustrated in figure Qa) . From figure |4jc) , we see that the error associated with ( 23 ) is of 
order C(n -1 / 2 ) in n. It should not come as a surprise that the leading term in (231 coincides with the spectrum of 
the Ornstein-Uhlenbeck process. However, it is not entirely obvious that the leading correction should also increase 
with n. Such a feature is explained by noting that the mismatch in width between the piecewise parabolic potential 
and the fully parabolic potential of the Ornstein-Uhlenbeck process increases in size when x increases. 
For large values of the parameter 8 and fixed values of n, the asymptotic property of the parabolic cylinder function 
results in 



AW + 1/2 



5 1 
4 



l+t„ 



-2/3 



2t 2 
15 



-4/3- 



«y> 1, 



(24) 



where — t n labels the zeros of the Airy function Ai(z). This expression, which is illustrated in figure |4]jb), captures the 
limiting case of vanishing viscous friction and the transition from a discrete to a continuous spectrum. One clearly 



FIG. 5: (Color online) Dry friction with viscous friction. Evolution of p(x, t\x' , 0) towards the stationary distribution (dashed 
line) as r increases. Parameters: x' = 2 and 5 = 1. 40 modes were used to plot the propagator for the two first times; for the 
the later times, only 5 modes were used. 



observes the spectral gap of size <5 2 /4, and the fact that the eigenvalues accumulate for large 5 in a quasicontinuous 
way when rescaled by 5 2 to account for the different time scales defined in ^ and (fl9l ). 

From figure [4^d), we are led to conjecture that the error term associated with (|24[) is of order 0{5~ 2 ). Note also 
that if we allow for negative S, then the asymptotic expansion of the parabolic cylinder function yields, as expected, 
an activation rate result for the Kramer escape problem, having the form 

A o) ~^i , (25) 



5 



with the denominator of the pre-exponential factor being given by the asymptotic series 

^g^(i)W-«e, f (|). 

Having studied the spectrum of the Fokker-Planck operator in detail, we now turn to the propagator. As in the 
pure dry friction case, the eigenfunctions can be classified as odd or even and are now given in terms of parabolic 
cylinder functions: 

u^l°\x)^e^{-{x + 5f/A)D^ lo) {x + 5), x>Q, (27) 

where the expression for negative x follows by symmetry. The eigenfunctions Vn^°^ (x) of the adjoint problem are, as 
usual, determined by The normalisation of the odd eigenfunctions reads 6 

4 o) = r v { : ] {x)u[f>{x)dx = 2k^D^_ x {5) Ud^D^iS)) , (28) 

J — OO 

and yields the normalisation of the even eigenfunctions via 

POO 

4 e) = / ^(x)u^(x) dx = (AW + 1)4°). (29) 



From these expressions, we obtain the spectral representation of the propagator as 

P (x,r\x',0) = Mx) + £ expl-AWr) ^^^ + £ exp^r) ""^ ( ^ . (30) 



This expression determines a nonnegative bounded real function. 

9aDa(S) stands for the derivative of the parabolic cylinder function with respect to its index. In addition to the standard derivation 
given in Appendix [B] the value of the normalisation integral can be computed by comparing the exact r esul t of the Laplace transform 
of the exit time distribution |27l with the Laplace transform of the odd part of the propagator shown in J30t. 




FIG. 6: (Color online) Dry friction with viscous friction, (a) Convergence of the mode decomposition of p(x,r\x' ,0) for 
n = 2,3,4 modes (coloured lines) and n = 20 modes (dashed black line). Parameters: x' = 3, 8 = 1 and r = 0.5. (b) Mode 
convergence for short time with n = 10, 50, 100, 200 modes. Parameters: x' = 3, 8 = 1 and r = 0.01. 



More details about this result can be found in Appendix [B] 

The behaviour of the propagator ( 30 1 for finite values of S is illustrated in figure [5j Overall, we see that for moderate 
values of the viscous damping, i.e., for moderate positive values of the parameter S, the evolution of p{x, t\x', 0) largely 
resembles the case of pure dry friction with its cusp at the origin; see figure jlj The only minor difference is that 
p(x,r\x' ,0) converges to the stationary distribution p*(x) relatively faster than with pure dry friction because of 
the added viscous friction. This property is reflected in the value of the lowest non-trivial eigenvalue at large 5; see 
equation p4| . 

In the absence of dry friction (S = 0), the potential $(x) becomes quadratic and the problem reduces to the 
Ornstein-Uhlenbeck process, as mentioned before. In this case, the parabolic cylinder functions for integer indices can 
be expressed in terms of Hermite polynomials 



D n (z) = e- z2 / 2 Re n (z), 



(31) 



and lead to odd and even eigenvalues corresponding to odd and even positive integers, respectively (see figure [3]). In 
this special case, the spectral decomposition (30) can be written in closed analytical form using Mehler's formula 25j. 



For the general case where dry and viscous frictions are present, we are not aware of a closed analytical expression 
for the propagator of (30). Thus, to evaluate this propagator, and to generate, for instance, the curves shown in figure 
[5j we have to rely on the numerical computation of the spectral sum. The behaviour of this sum as the number of 
eigenfunctions or modes is increased is shown in figure [6] As one expects from the result of ( 30 1 , one requires a large 
number of modes to capture the short time behaviour of p(x, t\x' , 0), and a Gibbs-type phenomenon is visible when 
an insufficient number of modes is used. This problem becomes irrelevant when the dynamics on larger time scales is 
considered, as modes in the upper part of the spectrum are strongly damped. 

For the correlation function, only the odd part of the propagator (30) contributes. Using ([9| for the adjoint 
eigenfunctions and ^ for the stationary distribution, we obtain the representation 



(x(t)x(O)) = 



y-r exp( 
n=l 



-Ai 0, r) 



HIo 



oo (o) 
XUn 



(x) dx 



Zk° ] JTexpH* + <5) 2 /2) dx 



(32) 



The numerator of this expression can be simplified using a recurrence relation for the parabolic cylinder functions 
reproduced in (B7) and (B8), the representation (27) for the eigenfunctions, and integration by parts to obtain 



x u 



(<>) 



(x) dx 



(o) 



AS). 



(33) 



Thus, the series (32) reads 



(x(t)x(O)) = 



4 exp(-(5 2 /2) ^ exp(-Al o) T) 
y/2n crfc((V\/2) 

S 2 



E 



7(0) 



-Di 



VV/2 exp((5 2 /2)erfc(<5/v / 2) 



E 



,(*) 



cx P (-A! o) t) 



a(°)(a(°) 

1 t^n \l±n 



l) 2 (" 



(34) 
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FIG. 7: (Color online) Dry friction with viscous fric tion. Log- linear plot of the normalised correlation function C(r) = 
{x(r)x(0))/{x 2 ) t for 5 



1, 2, 

oo; see (171 



obtained in the limit S 

case, the data are displayed in the units defined in |5| 



8, and 16; see |34j). The dashed line corresponds to the case without viscous damping, 
For the purpose of comparing the dry and viscous friction case with the pure dry friction 



where for the last step we have used (28) and (21), as well as the recurrence relation (B14). 

Figure [7] confirms that the correlation function ( 34 ) decays exponentially. This agrees with the spectral gap of the 
spectrum of the Fokker-Planck operator, which, according to equation (24), has a size of order <5 2 /4. For increasing 
values of 8, the correlation function of (34) approaches the pure dry friction result of de Gennes [T] shown in (17). In 
fact, for 5 ~ 10, the correlation function in the presence of dry and viscous frictions is almost indistinguishable from 
that obtained with dry friction only. From this, we conclude that a viscous damping of order 7r ~ 0.01A 2 does not 
affect de Gennes's result significantly. Since this result is based on the product 7r being small, we also conclude that 
de Gennes's result is a good approximation of the correlation function for the dry and viscous friction case when the 
noise is sufficiently weak. 

The effect of viscous damping becomes apparent only at large times, as a result of the discrete spectrum obtained 
with dry and viscous friction. The asymptotic result of (24) shows that the level spacing of the eigenvalues is of 
order S 2 ^ 3 , so that at a time scale of order r ~ <5~ 2 / 3 , we start to see a difference between the exponential decay of 
(x(t)x(O)) obtained with and without viscous damping. In the original units, this time scale is approximately given 
by 



i~i /(Wi ) 2/3 



(35) 



where i corr is the correlation time of the pure dry friction case, defined before in ( 18 1. The influence of viscous damping 
on the tail of the correlation function is thus visible only if "ft corr is not too large. One can study this regime in more 
detail by using asymptotic expansions for the parabolic cylinder function in the result of (34). 

To conclude this section, note that in the absence of dry friction, 6 = 0, the sum in (34) 
since the numerator D A ( a ) _ 1 (0) vanishes for n > 2 because of the characteristic equation 
expected, the correlation function of the Ornstein-Uhlenbeck process with viscous friction. 



contains only one mode, 
This mode recovers, as 



IV. CONSTANT EXTERNAL FORCE 



For a system with constant external force, the suitable nondimensional variables are again given by (19), and the 
potential &(x) which determines the drift of the Fokker-Planck equation reads 



(\x\ + Sf 



bx. 



A 



2 

0T 



1/2 



2 

tT 



1/2 



(36) 



The additional parameter b determines the strength of the "acceleration" a relative to the viscous damping. Depending 
on the relative size of the two parameters S and b, the deterministic part of the Langevin equation ([2 ) either corresponds 
to a particle at rest or to a moving particle, as mentioned in the introduction. For \b\ < 5, the deterministic steady 
state corresponding to the minimum of the potential (36) has vanishing velocity (stick state), while for \b\ > 6, a 
steady state with finite velocity occurs (slip state). In the presence of noise, these stick and slip steady states become 
the most probable stationary states of the Langevin dynamics, and so determine the maximum of the stationary 
distribution p*{x). This is illustrated in figure [8] 
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FIG. 8: (Color online) Dry friction with viscous friction and forcing, (a) Potential &(x) as a function of x for S 
different values of 6, as given by (36 1, (b) Corresponding invariant distribution, given by Q. 



1 and for 



To solve the time-dependent Fokker-Planck equation for the full-force case, we proceed along the lines of the previous 
sections. The main steps are given in Appendix [C] here we only summarise the main properties of the spectrum of 
the Fokker-Planck operator, the propagator p(x,r\x' ,0) and the corresponding correlation function. 

Unlike the previous case where 6 = 0, eigenfunctions cannot be classified according to their symmetry when 6^0. 
The characteristic equation is given by 



A n (D An (5 + b)D An ^(6 -b) + D An (5 - b)D An ^(6 + b)) = 0. 



(37) 



Without external forcing (b = 0), this expression reduces of course to the case discussed in the previous section; 
see (21) and (22). On the other hand, without dry friction (S = 0), we end up with the integer spectrum of the 
Ornstcin-Uhlenbeck process, since the Wronskian for the fundamental system (D„(z), D v (—z)) results in the identity 



D v {z)D v -x{-z) + D v (-z)D v -x(z) = V5F/T(-iy + 1). 



(38) 



Last but not least, we stress that the eigenvalues are even functions of the driving force b, since (37) is invariant under 
a change of sign. As a result, it is sufficient to consider the case of nonnegative driving b. 

Figure [9] illustrates the main properties of the spectrum for fixed positive values of 6. For large forcing, i.e., in 
the slip phase, the spectrum approaches the integer eigenvalues of the Ornstein-Uhlenbeck process, reflecting the 
parabolic minimum of the potential (see figure [8]). For small forcing and relatively small viscous damping, we observe 
the spectral gap between the ground state and the first nonzero eigenvalue, encountered before. In the transition 
regime between stick and slip phase, the eigenvalues show clear oscillations and a sharp drop when the transition 
value |6| = 5 is approached. The drop becomes more pronounced in the low noise limit, i.e., for large values of 5 and 6, 
and can be understood from (37) using our knowledge of parabolic cylinder functions gained in the previous section. 

To be more precise, consider the stick phase, i.e., < b < 8. As we have seen previously, the parabolic cylinder 
function D A (z) is positive if the index A is smaller than the lowest non-trivial eigenvalue of th e co rresponding 
undriven system, i.e., A < z 2 /4 for large positive argument z (see figure[2|. Thus all contributions to (37) are positive 
if A < (6 — 6) 2 /4, and this condition yields a lower bound for the lowest non-trivial eigenvalue, which determines, in 
turn, the parabolic shape of the spectrum in the stick phase visible in figure [9j This reasoning can be turned into an 
asymptotic expansion for the lowest eigenvalues if one recalls that a change in sign of one of the factors in ( 37 1 almost 
inevitably results in a zero of ( 37 1 , given that the amplitudes of the terms become exponentially large. 



In the slip phase b> 5, the potential (36) has a quadratic minimum (see figure [8]), and one expects the lower part 
of the spectrum to be of Ornstein-Uhlenbeck type. Such a feature can be derived from the characteristic equation 
(37). If the argument 5 — b in (|37| ) is negative, then the cylinder functions change sign at values for A close to even 
nonnegative integers (see, e.g., ( |25| for the lowest value or figure[3]). Since D A (6 + b) considerably exceeds D A _x(S + b) 
in magnitude, the first non-trivial root of (37) must then occur close to one, as observed in figure [9] By elaborating 
further on this argument, one can derive asymptotic expressions for the spectrum in the case 6 > b. 
We now turn to the eigenfunctions, which for nonvanishing eigenvalues are given by the expression 



i(x) = exp(-$0)/2) x 



XnD An {- 



b), 

-5 



b), 



x>0 
x<0 



(39) 
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FIG. 9: (Color online) Dry friction with viscous friction and forcing. First five eigenvalues A n as a function of b for 5 = 1 (a) 
and 5 — 5 (b). 




FIG. 10: (Color online) Dry friction with viscous friction and forcing. Evolution of p(x, t\x', 0) towards the stationary distri- 
bution (dashed line) as r increases for the case where b < 5. Parameters: x' = 3, b = 0.5, and 5=1. 



with the amplitude ratio 



D An (5-b) £>a„-i(5-&) 



(40) 



A " £> A „ (5 + 6) D A „-i(<5 + &)' 
The adjoint eigenfunctions v„(x) are as usual obtained from (j9j). The normalisation ( 11 ) can conveniently be expressed 



as 



Z n = A n D A JS - &)£>a„-i(<5 - b)d An In 



D Ab _ 1 (*-6)D Ab (5 + 6) 



£) a „(<5-6) j Da„-i(^ + &) 



Thus, (39 1, (40), and (41) provide the required input for the propagator 

u n {x)v n (x') 



p(x, t\x', 0) = p*{x) + ^2 exp(-A„r)- 



Z,, 



(41) 



(42) 



Figure 10 illustrates the time evolution of (42 ) in the stick phase |6| < S. Qualitatively, the dynamics is quite similar 
to the undriven system (see figure [5]), with the difference that the symmetry of the stationary state has disappeared. 

New features appear in the slip phase b > S > 0. Figure [TT] shows the time evolution of the propagator with two 
initial conditions of opposite sign. Depending on the sign of the initial condition relative to the stationary state, the 
propagator either remains unimodal or temporarily develops a bimodal shape. The latter is caused by a slow process 
near the origin, which was already responsible for the slow build-up of the tails of the stationary distribution in the 
undriven case. 
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FIG. 12: (Color online) Dry friction with viscous friction and 8 < 0. (a) First five eigenvalues as a function of b for S — — 1. 
(b) Evolution of p(x,r\x' ,0) towards the stationary distribution (dashed line) as r increases. Parameters: x — —3, S = — 1, 
and fe = 0.5. 



The bimodal form of p(x,r\x' ,0) is also seen for negative values of S, associated with Kramer's escape problem. 
Figure [12] gives an idea of the spectrum and the time evolution of the propagator in that case. One observes again 
a transition phenomenon in the spectrum when the shape of the potential changes from a bimodal to a unimodal 
structure. For 5 < 0, the driving force b affects the escape rate in the usual way. The lowest eigenvalue increases with 
increasing driving force as the potential barrier decreases in size. Moreover, the propagator displays, as expected, a 
slow tunnelling process in a double-well potential. 

A few more insights into the forced case can be obtained by computing the correlation function along the lines of 
the previous section; see, e.g., (32) and (33). In the present case, (39) yields, using (B7) and integration by parts, 



xu n (x)dx = e- s2 / 4 [D A ^ 2 (S-b)-XAD An - 2 (S + b)} 



A, 



1 



-DA^iS-b), 



(43) 



where for the last step the recurrence relation (B14) and the definition (40) have been used. From the propagator 
(42), we then obtain 



(x(r)x(O)) - (x)l - 



AS 2 



,(«-6)V2 erfc (*jf) + e (W/2 erfc 
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FIG. 13: (Color online) Dry friction with viscous friction and forcing. Normalised correlation function C(r) = ((x(t)x(O)) — 
{x)1) / ({x 2 ) * — (x)t) as a function of r and b for S — 12. The red line at 6 = 5 marks the transition between the stick and slip 
states. 



exp( 



-A* 



y 

£j K(An-l) 2 D An (8-b)d An ln 



D An ^(d-b) 



D An -i(S-b)D An (S+b) 
£>A„(<S-6)£>A„-i(5+b) 



(44) 



In the limit b — > 0, this result reduces to the correlation function of the undriven case, (34), as the denominator of the 



expansion coefficient tends towards — 2d\ n D An (8) and the numerator ensures that only the odd part of the spectrum 
contributes; see (22). 



Figure [13] shows the plot of the correlation function normalised by the stationary variance as a function of r and 
b. The red line in this plot marks the boundary b = 5 separating the noiseless or deterministic stick and slip states. 
In the presence of noise, these two states are now recognisable as two different correlation phases. In the stick phase, 
< b < 5, it can be shown from the spectrum of the Fokker-Planck operator that the correlation time governing 
the exponential decay of (x(t)x(O)) scales according to (8 — b)~ 2 , and so dramatically increases near the stick-slip 
transition. In the slip phase (b > 8), by contrast, the correlation time is mostly determined by the viscous damping, 
and thus does not change sensibly as a function of b. This difference in the behaviour of (x(t)x(O)) is clearly seen in 
figure [F3| and becomes more visible in the small- noise limit. 



V. CONCLUSIONS 



In this paper we have studied a Langevin equation that includes a solid or dry friction force and an external constant 
force in addition to the viscous friction force commonly considered in the context of Brownian motion. This stochastic 
equation extends an earlier model studied by de Gennes [1] , which can be thought of as describing the dynamics of a 
solid object resting on a tilted solid surface, which is vibrated randomly with Gaussian white noise. By solving the 
time-dependent Fokker-Planck equation associated with this model, we have obtained the time-dependent propagator, 
which gives the probability that the object has a certain velocity at any time starting from a given initial velocity, in 
addition to the velocity correlation function. 

These results, combined with the full spectrum of the Fokker-Planck operator, allowed us to obtain a detailed 
understanding of dry friction in the presence of noise. In particular, we have seen that the singular nature of the 
dry friction force at zero velocity gives rise to a cusp in the propagator whenever the external forcing is smaller in 
magnitude than the dry friction force. We have also seen that the viscous damping does not alter the properties of 
the model much when the noise is sufficiently small, essentially because the stochastic dynamics is then confined near 
the origin (zero velocity state) which is most sensitive to dry friction. Conversely, for large noise powers, the model is 
mostly dominated by diffusion and behaves similarly as an Ornstein-Uhlenbeck process. Finally, we have seen that, 
although the stick-slip transition is blurred at the level of the stochastic dynamics, it is possible to define stick and slip 
phases at the level of the correlation function. The stick phase is characterized by a sharp increase of the correlation 
time with the external force, whereas, for the slip phase, the correlation time is mostly determined by diffusion and 
is almost independent of the external force. These quantitative predictions could be checked in real experiments and 
in more elaborate models of dry friction. 
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To conclude, it is worth noting that the results presented here have direct analogs in terms of quantum mechanical 
potentials. Indeed, some of our results about the spectrum of the Fokker-Planck operator in the presence of dry 
and viscous frictions can be inferred from previously-published results on the so-called constrained quantum harmonic 
oscillator [28 , 29 . However, to the best of our knowledge, the solution that we have obtained for the full- force case has 
not been considered before. It might be of interest to translate this solution to the quantum context. Furthermore, 
it might be interesting, at least from an aesthetic point of view, to see whether the spectral sums obtained for dry 
and viscous forces and the full-force case can be written in closed analytical forms, as in the case of pure dry friction. 
This might be possible given that the Fokker-Planck equation is piecewise linear in all cases. 
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Appendix A: Particle in a wedge 



The stochastic dynamics of a particle in a symmetric wedge potential 



(Al) 



is a classical textbook example of a Fokker-Planck operator with spectral gap and continuous spectrum [24]. Here 
we recall the solution of this well-known and somewhat elementary problem, as it will certainly help the interested 
reader to understand the more elaborate considerations found in the other appendices. 



Eigenfunctions are either even, u A (— x) — u^'{x), or odd, 



(<0, 



-x) 



-u A ' (x), with respect to inversion. The 



eigenvalue equation ([8|) may be considered on the positive domain only and the boundary condition, given before in 
(13), readily reduces to 



u w (0)J *4 e) (*) 



dx 



< ) (0)=0. 



(A2) 



The general solution of (JsJ) can be expressed as a linear combination of exponentials. Apart from the stationary 
distribution ^ with eigenvalue A = 0, the solution which satisfies the symmetry, the boundary condition (A2), and 
which decays to zero at infinity is given by 



«i e) 0O 
u<°\x) 



e~ x ' 2 [2/3 A cos(/3 A :r)-sin(/3 Aa; )] 
e~ x ' 2 sm(j3 A x) 



for x > 0, where 



(3 A = - 1/4 > 0. 



(A3) 



(A4) 



The latter inequality determines the spectral gap and restricts the continuous spectrum to the range A > 1/4. As for 
the normalisation of the odd eigenmodes, (A3) and Q result in 



v A )(x) u A (x) dx 



whereas for the even eigenmodes, we obtain 



[cos ((/9a - Pa>)x) - cos((/3 A + /3a>)x)} dx 
2ny/A-l/4S(A - A') = Z ( A o) 8(A - A'), 



(A5) 



v a? (x) uf^ (x) dx 



[(4/? A /9a' + 1)cos(C9a-Av)z) 

+(4/?a/3a< - 1) cos((/3 A + fa)x) + 2(/9 A - A ,) sm((/3 A - /3 A ,)x) 

-2(/3 A + p K ,) sin((/9 A + Mx)] dx 
2ir(4(3 2 A + l)v/A- 1/4<J(A - A') 
Z A e) <5(A- A'). 



(A6) 
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The contribution of the odd eigenmodes to the propagator (10) can now be evaluated straightforwardly using (A3 1, 
@ and ( |A5| ): 



P (°\x,t\x',0) = 



e- AT u ( °\x)v { °\x')/Z { ° ) dA 



1/4 

e -r/i e -(\x\-\x'\)/2 r v 
27T 

e -r/i e -(\x\-\x'\)/2 



e-' 3 T [cos{f3{x - a;')) - cos(/3(:r + x'))} d/3 



e -(x-x') 2 /(4r) _ e -(x+x') 2 /(±r) 



For the even modes, we obtain instead 
p (e) (^,r|a;',0) = 



1/4 



ut\x)vt\x')/Z^dk 



r/± p -(\x\-\x'\)/2 /-oo -/3 



2tt 



+ (4/3 2 - 1) cos(p{\x\ + \x'\)) - 4/3sin(/3(|x| + \x'\))] d(3. 
The remaining non-trivial integrals can be evaluated using the identities |34j 



-/3x 2 



sin(ax) 



~P X " cos{ax) 



x dx 




7 2 + x 2 


= ~v 


dx 


7T 


7 2 + x 2 


47 



Pr 



2cosh(a 7 ) - e _7a erf ( 27/3 _° ) - < 'crl' 



We thus end up with 7 

e -T/i e -{\x\-\x'\)/2 



p {e \x,T\x',0) 



4\/7rr 



; -(x-o;') 2 /(4r) + e -(x+x') 2 / '(4r) 



2V? 



-kl 



2V? 



1 — erf 



2 V* 



The final result, ( 14 1, then follows from (A7), (All) and 



p(x, t\x', 0) = p, (x) + V [e) (x, t\x', 0) + p (o) (sb, t\x', 0). 
Appendix B: Symmetric piecewise parabolic potential 



(A7) 



(A8) 

(A9) 
(A10) 

(All) 
(A12) 



For the symmetric potential ( 20 ) , it is sufficient to consider the eigenvalue equation ([8|) for nonnegative argument 
with appropriate boundary condition at the origin, since eigenfunctions are either even or odd. The boundary condition 
expressed in (13) here yields (see (A2|) 



<H (o) 



dx 



= 0, z4 o) (0) = 0. 



(Bl) 



x=Q 



The differential equation (|8| with the potential (20) is a special case of the so-called Kummer's equation, whose 
solution can be expressed in terms of Kummer functions. The solution which decays at infinity can conveniently be 
expressed in terms the parabolic cylinder functions. 8 To be more precise, we have that 



solves 



u{z) = e- z2 / 4 D„(z), 
u" + (zu)' + vu = 0, 



(B2) 
(B3) 



7 See [26 for a closely related result. 

8 We use Whittaker's notation; see, e.g., |30| 
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and D v [z) ~ e z / 4 as z — > +00. Thus the solution of (|Sp , which decays at infinity, is given by 

u A (x) =exp(-(.x + (5) 2 /4)L>A(a; + '5), a; > 0. (B4) 

The spectrum is now determined by the boundary condition at the origin. For the odd eigenf unctions, ( |B1[ ) and 
(B4| yield the condition (21 ). If we use the following differential identity for the parabolic cylinder function: 



D' A (z) + -D A (z) = AD A _ 1 (z), 
then the boundary condition for the even eigenfunctions results in 



(B5) 



(B6) 



This means that, apart from the trivial zero eigenvalue associated with the stationary state, all the other eigenvalues 
obey (22) if we take the result of (21 ) into account. 

Turning now to the eigenfunctions, it should be clear that (B4| leads to (27). The adjoint eigenfunctions are 
determined by (J9|. The orthogonality of the eigenfunctions follows from the standard differential identities for the 
parabolic cylinder function, which are expressed in (|B5[) and can be written as 



<iexp(-z 2 /4)Z)„(z) 
dz 

de*v(z 2 /4)D v {z) 
dz 



= -cxp(-z 2 /4)A, +1 (z) 
= z/exp(z 2 / 4 )A>-i(»- 



As for the normalisation, using (B7) and <|B8j) , we obtain the relations 

dv ^\ X ) _ _„(e) (x) dV ^( X ) _ A (e) v (o) (x) 

dx ~ n 1 >' dx ~ n n [ h 



(B7) 
(B8) 

(B9) 



relating the eigenfunctions and the eigenfunctions of the adjoint problem. Integration by parts thus leads to 

(x) «#> (x) dx = / »W (x) U W ( X ) dx. (BIO) 



This last equation yields (29). In fact, (B9) are the remnants of the differential relations for Hermite polynomials, 
which guarantee that for the case 5 = 0, i.e., for the Ornstein-Uhlenbeck process, all eigenfunctions can be obtained 
by the shift operation. 

To evaluate the normalisation integrals and to establish the orthogonality of eigenfunctions, we resort to a known 
integral identity of parabolic cylinder functions which follows from (B7) and (B8) using integration by parts [33] : 



(A -A') / D A (x)D A ,(x)dx 



dexp(x 2 /4)£>A+i(» 



dx 



exp(a; 2 /4)D A (x) 



e*v{-x 2 /A)D A ,{x) 
dexp(x 2 /4)D A , +1 (x) 



dx 



dx 



= D A (a)D A , +1 (a) - D A+1 (a)D A ,(a). 



(Bll) 



Choosing a = 5 and A^ ' = A 7^ A' = A^', the right hand side of (Bll) vanishes due to the characteristic equation 
(21 ) and the left hand side results in the orthogonality condition 

/>oo 

(x) i4 o) (x) dx = 2 D aM (x + 6)D A(o) (x + S)dx = 0, 



» 



m =/= n 



(B12) 



for parabolic cylinder functions with index determined by (21). This relation is expected, as eigenfunctions are 
supposed to be mutually orthogonal. 



To evaluate the normalisation integral, let us consider the case a = S and A = A„ \ Then (Bll ) results in 



D 



A (o)(x)D A ,(x)dx = D A (o) +1 (5) 



D Aln0) (S)-DA>(S) 

A' 



(B13) 



if we again take the characteristic equation into account. In the limit A' — > An \ we then obtain the result (28) using 
the recurrence relation 



zD v {z) - D u+1 {z) - vD v - X {z) = . 



(B14) 
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Appendix C: Biased forcing 

For the potential (36 1, the eigenvalue equation Q reads 

u" A (x) + ((x±S- b)u A )' + Au A (x) = 0, 



(CI) 



where the + sign applies for x > and the — sign applies for x < 0. The solution that fulfills the boundary condition 
at infinity can be written again in terms of parabolic cylinder functions: 



u A ( x ) = C± exp(-(x ±8- b) 2 /A)D A (±x + STb), x^O. 



The boundary conditions ( 13 1 thus result in 

C_ exp(-(<5 + b) 2 /4)£> A (5 + b) 
C_Aexp(-(<5 + 6) 2 /4)D A _ 1 (5 + 6) 



C + exp(-(S-b) 2 /4)D A (8~b) 
-C + Aexp(-(<5 - b) 2 /A)D A _ 1 (8 - b). 



(C2) 



(C3) 
(C4) 



The condition for a non-trivial solution yields the characteristic equation (|37j). As for the amplitudes C±, there 
does not seem to be a convenient form which can cope both with the limit of sym m etric poten tials and the diff erent 
symmetries of the eigenf unctions. If we choose C+ = exp(6 2 /4 — (56/2), then (C2|, (C3), and (C4| lead to (39) and 

@- i-i i-i i-i 

The normalisation (11) of the eigenfunctions (39) and 191) reads 



Z n — 



D An (x + S- b)D An (x + 8- b)dx + \ n \ D A„ (x + 8 + b)D An (x + 8 + b)dx. 



(C5) 



Using the identity (Bll ) with a = 8 ± b, A = A„ and A' — > A„, we can express the remaining integrals in terms of 
parabolic cylinder functions to obtain 



D An+1 (S - b)d An D An (S -b)- D An (5 - b)d An D An+l {8 - b) 

+X 2 n [D An+1 {6 + b)d An D An {8 + b) - D An {8 + b)d An D An+1 (8 + b)} 



(C6) 



Using the identity (B14) to eliminate the cylinder functions with the largest index, D An+ i(S ± 6), and applying the 
definition (40), we finally arrive at the result of (41). 
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